This project investigates the complex relationship between transit ridership and public safety around SEPTA bus stops in Philadelphia. By employing Negative Binomial Regression, we explore the tension between two competing theories: whether high passenger volumes act as protective “Eyes on the Street” or attract crime as “Potential Targets.” The analysis integrates multi-source data—including ridership, crime incidents, census demographics, and built environment features—and utilizes interaction terms to reveal how risk dynamics shift between weekdays and weekends. Our ultimate goal is to transition SEPTA Transit Police from reactive to predictive deployment, providing a data-driven framework to identify high-risk anomalies and optimize patrol resources for maximum safety impact.
The primary dataset we use is 2025summer Septa bus ridership data.The validity of any predictive model rests on the quality of its underlying data. In this initial phase, we focus strictly on data preparation and exploratory visualization. We load and clean the raw ridership, crime dataset and other datasets to handle missing values and outliers. Through initial descriptive analytics and visualization, we assess the statistical properties of our datasets.
Load necessary libraries
library(tidyverse)
library(sf)
library(tidycensus)
library(tigris)
options(tigris_use_cache = TRUE, tigris_class = "sf")
library(MASS)
library(spdep)
library(dplyr)
library(scales)
library(ggplot2)
library(caret)
library(nngeo)
library(car)
library(knitr)
library(readr)
library(patchwork)
library(kableExtra)
library(lubridate)
options(tigris_use_cache = TRUE)
options(tigris_progress = FALSE) Define themes
plotTheme <- theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(size = 8),
axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
axis.title = element_text(size = 11, face = "bold"),
panel.background = element_blank(),
panel.grid.major = element_line(colour = "#D0D0D0", size = 0.2),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
legend.position = "right"
)
mapTheme <- theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(size = 8),
axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(colour = 'transparent'),
panel.grid.minor = element_blank(),
legend.position = "right",
plot.margin = margin(1, 1, 1, 1, 'cm'),
legend.key.height = unit(1, "cm"),
legend.key.width = unit(0.2, "cm")
)
palette5 <- c("#eff3ff", "#bdd7e7", "#6baed6", "#3182bd", "#08519c")2025 summer SEPTA bus ridership data
# 1. Load Bus Data
bus_raw <- read_csv("data/Summer_2025_Stop_Summary_(Bus).csv")# 2. Get Philadelphia Boundary
philly_boundary <- counties(state = "PA", cb = TRUE, class = "sf") %>%
filter(NAME == "Philadelphia") %>%
st_transform(2272)In preparing the ridership data, we performed two critical structural transformations to align with our research goals.
First, we aggregated bidirectional stops sharing the same location (including opposite sides of a street) into single spatial units. This prevents spatial autocorrelation and ensures our 400m buffers capture a unique street environment.
Second, we restructured the dataset into a ‘long format’ (panel data), distinguishing between Weekday and Weekend ridership. This temporal split is vital, as it allows our model to capture the distinct behavioral patterns of both commuters and potential offenders during different times of the week.
# 3. Process Ridership: Create Long Format (Aggregated + Weekday vs Weekend)
bus_long <- bus_raw %>%
filter(!is.na(Lat) & !is.na(Lon)) %>%
mutate(
Raw_Weekday = Weekdays_O + Weekdays_1,
# A. Weekend average = (Sat total + Sun total) / 2
Raw_Weekend = (Saturdays_ + Saturdays1 + Sundays_On + Sundays_Of) / 2
) %>%
# B. Aggregation by Stop Name
group_by(Stop) %>%
summarise(
# C. Add two direction ridership
Ridership_Weekday = sum(Raw_Weekday, na.rm = TRUE),
Ridership_Weekend = sum(Raw_Weekend, na.rm = TRUE),
Lat = mean(Lat, na.rm = TRUE),
Lon = mean(Lon, na.rm = TRUE),
# Keep one Stop_Code as unique ID
Stop_Code = first(Stop_Code),
.groups = "drop"
) %>%
# D. Pivot to Long Format
pivot_longer(
cols = c(Ridership_Weekday, Ridership_Weekend),
names_to = "Time_Type",
values_to = "Ridership"
) %>%
# E. Create Dummy Variable & Clean up
mutate(
is_weekend = if_else(Time_Type == "Ridership_Weekend", 1, 0),
Time_Category = if_else(is_weekend == 1, "Weekend", "Weekday")
) %>%
# F. Convert to SF & Clip
st_as_sf(coords = c("Lon", "Lat"), crs = 4326) %>%
st_transform(2272) %>%
st_intersection(philly_boundary) %>%
dplyr::select(Stop_Code, Stop, Ridership, is_weekend, Time_Category, geometry)
# Check results
cat("Total Aggregated Stops (Unique Locations):", nrow(bus_long) / 2, "\n")Total Aggregated Stops (Unique Locations): 5884 cat("Total Observation Rows (Panel Data):", nrow(bus_long), "\n")Total Observation Rows (Panel Data): 11768 ggplot() +
# A. Base Layer: Philadelphia County Background
geom_sf(data = philly_boundary,
fill = "grey98",
color = "grey50",
size = 0.5) +
# B. Station Layer: Simple Points
geom_sf(data = bus_long %>% filter(is_weekend == 0),
color = "#3182bd", # SEPTA Blue
size = 0.1, # Small dots to avoid clutter
alpha = 0.6) + # Slight transparency
# C. Styling
labs(
title = "Spatial Coverage of SEPTA Bus Network",
subtitle = paste0("Total Aggregated Stops: ", nrow(bus_long) / 2),
caption = "Source: SEPTA Summer 2025 Stop Summary"
) +
theme_void() + # Clean look
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 10, color = "grey40", hjust = 0.5),
plot.margin = margin(1, 1, 1, 1, "cm")
)
The figure illustrates the extensive spatial coverage of the SEPTA bus network within Philadelphia. After aggregating bidirectional stops, our study includes 5,884 unique locations.
The distribution reveals a high density in the Center City and a grid-like arterial pattern extending into residential neighborhoods. This comprehensive coverage ensures that our analysis captures a diverse range of built environments, from dense commercial districts to suburban residential areas.
We excluded domestic disputes or indoor financial crimes because SEPTA Transit Police cannot police inside people’s homes. Instead, we selected these specific crime types because they occur in the public right-of-way and directly impact a rider’s decision to use public transit.
### 1.2.1 Crime data:
# Define Street-Crime List
# target_crime_types <- c(
# "Robbery Firearm",
# "Robbery No Firearm",
# "Aggravated Assault Firearm",
# "Aggravated Assault No Firearm",
# "Homicide - Criminal",
# "Rape",
# "Other Sex Offenses (Not Commercialized)",
# "Weapon Violations",
# "Narcotic / Drug Law Violations",
# "Vandalism/Criminal Mischief",
# "Prostitution and Commercialized Vice",
# "Public Drunkenness"
# )
crime_raw <- read_csv("data/crime2024.csv")
# Filter & Transform
crime_raw <- crime_raw %>%
filter(!is.na(lat) & !is.na(lng))
crime_q1 <- crime_raw %>%
filter(lubridate::month(dispatch_date) %in% 1:3) %>%
mutate(
crime_date = as.Date(dispatch_date),
day_of_week = wday(crime_date), # 1 is Sunday, 7 is Saturday
is_crime_weekend = if_else(day_of_week %in% c(1, 7), 1, 0)
) %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_transform(2272)
# Filter Apr to Jun
crime_sf <- crime_raw %>%
filter(lubridate::month(dispatch_date) %in% 4:6)%>%
#filter(text_general_code %in% target_crime_types) %>%
mutate(
crime_date = as.Date(dispatch_date),
day_of_week = wday(crime_date), # 1 is Sunday, 7 is Saturday
is_crime_weekend = if_else(day_of_week %in% c(1, 7), 1, 0)
) %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_transform(2272)
cat("Total Selected Crimes (Count):", nrow(crime_sf), "\n")Total Selected Crimes (Count): 24184 To isolate the true impact of transit ridership on crime, we must first account for the socio-economic context of the neighborhood. A bus stop in a distressed area faces different risks than one in a wealthy suburb.
By integrating 2023 ACS Census data, we control for structural disadvantages—such as poverty rates, unemployment, and housing vacancy.
census_api_key("42bf8a20a3df1def380f330cf7edad0dd5842ce6", overwrite = TRUE, install = TRUE)# Load Census data for Philadelphia tracts
philly_census <- get_acs(
geography = "tract",
variables = c(
total_pop = "B01003_001",
poverty_pop = "B17001_002",
med_income = "B19013_001",
ba_degree = "B15003_022",
total_edu = "B15003_001",
labor_force = "B23025_003",
unemployed = "B23025_005",
total_housing = "B25002_001",
vacant_housing = "B25002_003"
),
year = 2023,
state = "PA",
county = "Philadelphia",
geometry = TRUE,
output = "wide"
) %>%
st_transform(2272) %>%
mutate(
Poverty_Rate = poverty_popE / total_popE,
Med_Income = med_incomeE,
ba_rate = 100 * ba_degreeE / total_eduE,
unemployment_rate = 100 * unemployedE / labor_forceE,
vacancy_rate = 100 * vacant_housingE / total_housingE
) %>%
dplyr::select(GEOID, Poverty_Rate, Med_Income, ba_rate, unemployment_rate, vacancy_rate)We ingest data from OpenDataPhilly to capture the built environment’s impact on safety. This includes Alcohol Outlets (potential crime generators), Street Lights (natural surveillance/visibility) and Police Stations (formal guardianship/deterrence).
# A. Alcohol Outlets (Crime Generators)
alcohol_sf <- read_csv("data/PHL_PLCB_geocoded.csv") %>%
filter(!is.na(lon) & !is.na(lat)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(2272)# B. Street Lights (Guardianship)
lights_sf <- read_csv("data/Street_Poles.csv") %>%
filter(!is.na(X) & !is.na(Y)) %>%
st_as_sf(coords = c("X", "Y"), crs = 3857) %>%
st_transform(2272)Proximity to police stations acts as a proxy for formal guardianship. We load these coordinates to later calculate the distance-based deterrence effect.
# C. Police Stations (Guardianship)
police_sf <- st_read("data/Police_Stations.geojson", quiet = TRUE) %>%
st_transform(2272)
cat("Total Police Stations:", nrow(police_sf))Total Police Stations: 20Phase 2 transforms our raw spatial data into analytical features. A bus stop is not just a point on a map; it is the center of a dynamic micro-environment. To capture this, we employ a 400-meter buffer approach—roughly equivalent to a 5-minute walk, defining the ‘catchment area’ for each station.
After that, we shift our focus to investigating the aggregate relationships between our independent variables (e.g., ridership, built environment) and the dependent variable (crime counts). This exploratory analysis allows us to observe baseline correlations and structural patterns, providing the necessary empirical support for our subsequent Negative Binomial modeling.
# Create Buffer(400m)
bus_buffer <- st_buffer(bus_long, 1312)# 1. Calculate total count
# Calculate crime_sf data set covers how many weekdays and weekends
day_counts <- crime_sf %>%
st_drop_geometry() %>%
group_by(is_crime_weekend) %>%
summarise(n_days = n_distinct(crime_date))
print(day_counts)# A tibble: 2 × 2
is_crime_weekend n_days
<dbl> <int>
1 0 65
2 1 26# 2. Join and calculate daily crime Rate
crime_agg <- st_join(bus_buffer, crime_sf, join = st_intersects) %>%
filter(is_weekend == is_crime_weekend) %>%
group_by(Stop_Code, is_weekend) %>%
summarise(
Crime_Total_Count = n()
) %>%
st_drop_geometry() %>%
left_join(day_counts, by = c("is_weekend" = "is_crime_weekend")) %>%
mutate(
Crime_Daily_Rate = Crime_Total_Count / n_days,
# Keep n_days，as negative binomial regresssion offset
Exposure_Days = n_days
)
head(crime_agg)# A tibble: 6 × 6
# Groups: Stop_Code [3]
Stop_Code is_weekend Crime_Total_Count n_days Crime_Daily_Rate Exposure_Days
<dbl> <dbl> <int> <int> <dbl> <int>
1 2 0 11 65 0.169 65
2 2 1 2 26 0.0769 26
3 4 0 70 65 1.08 65
4 4 1 26 26 1 26
5 5 0 102 65 1.57 65
6 5 1 34 26 1.31 26We know that a dangerous corner doesn’t become safe overnight. To make our model fair, we must acknowledge that some stops have a history of trouble.
We calculated the crime counts for each stop from the previous season (Q1). We try to explain: ‘Given how dangerous this spot usually is, does adding more bus riders make it safer or worse?’
day_counts_q1 <- crime_q1 %>%
st_drop_geometry() %>%
group_by(is_crime_weekend) %>%
summarise(n_days_lag = n_distinct(crime_date))
print(day_counts_q1)# A tibble: 2 × 2
is_crime_weekend n_days_lag
<dbl> <int>
1 0 65
2 1 26crime_lag_q1 <- st_join(bus_buffer, crime_q1, join = st_intersects) %>%
filter(is_weekend == is_crime_weekend) %>%
group_by(Stop_Code, is_weekend) %>%
summarise(
Crime_Total_Count_lag = n()
) %>%
st_drop_geometry() %>%
left_join(day_counts_q1, by = c("is_weekend" = "is_crime_weekend")) %>%
mutate(
Crime_Daily_Rate_lag = Crime_Total_Count_lag / n_days_lag,
Exposure_Days_lag = n_days_lag
)
head(crime_lag_q1)# A tibble: 6 × 6
# Groups: Stop_Code [3]
Stop_Code is_weekend Crime_Total_Count_lag n_days_lag Crime_Daily_Rate_lag
<dbl> <dbl> <int> <int> <dbl>
1 2 0 5 65 0.0769
2 2 1 2 26 0.0769
3 4 0 63 65 0.969
4 4 1 26 26 1
5 5 0 52 65 0.8
6 5 1 27 26 1.04
# ℹ 1 more variable: Exposure_Days_lag <int>Alcohol outlets are established ‘crime generators.’ We load their locations to model areas with higher potential for intoxication-related conflicts.
alcohol_agg <- st_join(bus_buffer, alcohol_sf, join = st_intersects) %>%
group_by(Stop_Code, is_weekend) %>%
summarise(Alcohol_Count = n() - 1) %>% # Subtract 1 because st_join is left join (self-intersection NA check)
st_drop_geometry()Street lighting is a key component of CPTED (Crime Prevention Through Environmental Design). We process the location of street poles to estimate visibility and natural surveillance levels.
light_agg <- st_join(bus_buffer, lights_sf, join = st_intersects) %>%
group_by(Stop_Code, is_weekend) %>%
summarise(Light_Count = n() - 1) %>%
st_drop_geometry()To control for neighborhood context, we calculated the average:
- Avg_Poverty: Rate of economic deprivation.
- Avg_Income: Proxy for local wealth and resources.
- Avg_BA: Educational attainment.
- Avg_Unemployment: Measure of labor market instability.
- Avg_Vacancy: Indicator of physical disorder.
# Average Census Demographics
census_agg <- st_join(bus_buffer, philly_census, join = st_intersects) %>%
group_by(Stop_Code, is_weekend) %>%
summarise(
Avg_Poverty = mean(Poverty_Rate, na.rm = TRUE),
Avg_Income = mean(Med_Income, na.rm = TRUE),
Avg_BA = mean(ba_rate, na.rm = TRUE),
Avg_Unemployment = mean(unemployment_rate, na.rm = TRUE),
Avg_Vacancy = mean(vacancy_rate, na.rm = TRUE)
) %>%
st_drop_geometry()# Calculate distance to nearest police station
dist_matrix <- st_distance(bus_long, police_sf)
bus_long$Dist_Police <- apply(dist_matrix, 1, min) / 5280psa_sf <- st_read("data/Boundaries_PSA.geojson", quiet = TRUE) %>%
st_transform(2272) %>%
dplyr::select(PSA_ID = PSA_NUM)
# Spatial Match: Define each bus stop to PSA
stop_psa_mapping <- bus_long %>%
filter(is_weekend == 0) %>%
st_join(psa_sf) %>% # Spatial Join：Point to Polygon
st_drop_geometry() %>%
dplyr::select(Stop_Code, PSA_ID) %>%
distinct(Stop_Code, .keep_all = TRUE)
cat("PSA mapping created for", nrow(stop_psa_mapping), "stops.\n")PSA mapping created for 5884 stops.head(stop_psa_mapping)# A tibble: 6 × 2
Stop_Code PSA_ID
<dbl> <chr>
1 759 033
2 16397 251
3 17921 351
4 16054 352
5 16102 221
6 16104 221 final_data <- bus_long %>%
# Join all aggregated tables
left_join(crime_agg, by = c("Stop_Code", "is_weekend")) %>%
left_join(crime_lag_q1, by = c("Stop_Code", "is_weekend")) %>%
left_join(alcohol_agg, by = c("Stop_Code", "is_weekend")) %>%
left_join(light_agg, by = c("Stop_Code", "is_weekend")) %>%
left_join(census_agg, by = c("Stop_Code", "is_weekend")) %>%
left_join(stop_psa_mapping, by = "Stop_Code") %>%
mutate(
# 1. Handle NAs for Counts
Crime_Total_Count = replace_na(Crime_Total_Count, 0),
Crime_Daily_Rate = replace_na(Crime_Daily_Rate, 0),
# Deal with NA in lag crime
Crime_Total_Count_lag = replace_na(Crime_Total_Count_lag, 0),
Crime_Daily_Rate_lag = replace_na(Crime_Daily_Rate_lag, 0),
Alcohol_Count = replace_na(Alcohol_Count, 0),
Light_Count = replace_na(Light_Count, 0),
# 2. Log transform Ridership
Log_Ridership = log(Ridership + 1),
# 3. Create Factor for Interaction Term
is_weekend_factor = factor(is_weekend, levels = c(0, 1),
labels = c("Weekday", "Weekend")),
) %>%
# 4. Clean up
filter(!is.na(PSA_ID)) %>%
na.omit()
cat("Final Panel Dataset Rows:", nrow(final_data))Final Panel Dataset Rows: 11066# 1. Panel Data Rows
psa_counts <- final_data %>%
group_by(PSA_ID) %>%
summarise(
n_observations = n(),
n_stops = n_distinct(Stop_Code)
) %>%
arrange(n_observations)
print(head(psa_counts, 10))Simple feature collection with 10 features and 3 fields
Geometry type: MULTIPOINT
Dimension: XY
Bounding box: xmin: 2669882 ymin: 215531.6 xmax: 2715062 ymax: 279357.4
Projected CRS: NAD83 / Pennsylvania South (ftUS)
# A tibble: 10 × 4
PSA_ID n_observations n_stops geometry
<chr> <int> <int> <MULTIPOINT [US_survey_foot]>
1 122 32 16 ((2670537 223892.4), (2670585 223932.8), (2670…
2 124 58 29 ((2671643 232878.2), (2672084 233106.1), (2672…
3 053 68 42 ((2670063 279330.9), (2670112 279357.4), (2670…
4 012 83 44 ((2685082 218564.9), (2685568 218997.3), (2686…
5 224 86 43 ((2686092 246848.8), (2686176 246409.7), (2686…
6 393 90 45 ((2688934 252812), (2688968 253070.2), (268922…
7 181 96 48 ((2669882 237256), (2669902 237214.7), (266998…
8 243 96 49 ((2706035 247307.4), (2706350 248302.6), (2706…
9 242 98 49 ((2702265 250424.1), (2702297 250520.3), (2702…
10 052 99 50 ((2671211 266842), (2671303 266812.6), (267140…Does ridership follow a normal distribution? No. That’s why we need Negative Binomial.
library(gridExtra)
p1 <- ggplot(final_data, aes(x = Ridership)) +
geom_histogram(fill = "#3182bd", bins = 50) +
labs(title = "Distribution of Ridership", x = "Daily Boardings") + plotTheme
p2 <- ggplot(final_data, aes(x = Crime_Total_Count)) +
geom_histogram(fill = "#de2d26", bins = 50) +
labs(title = "Distribution of Crime", x = "Crime Count (400m)") + plotTheme
grid.arrange(p1, p2, ncol = 2)
The histograms reveal that both Ridership (Left) and Crime Counts (Right) follow a highly right-skewed, ‘long-tail’ distribution, rather than a normal bell curve. This extreme skewness mathematically confirm that a standard OLS Linear Regression would be biased.
This visual evidence creates a compelling mandate for using Negative Binomial Regression, which is specifically designed to handle such skewed count data.
mean_crime <- mean(final_data$Crime_Total_Count, na.rm = TRUE)
var_crime <- var(final_data$Crime_Total_Count, na.rm = TRUE)
cat("Mean:", mean_crime, "\n")Mean: 30.91478 cat("Variance:", var_crime, "\n")Variance: 1190.485 cat("Ratio (Var/Mean):", var_crime / mean_crime, "\n")Ratio (Var/Mean): 38.50861 While many bus stops have zero incidents, a few ‘hotspots’ have very high counts. This causes the variance to be much larger than the mean. To address this overdispersion，which violates the core assumption of Poisson regression， we employed a Negative Binomial model. This approach allows us to model the data more accurately without inflating the significance of our findings.
# 1. Extract Coordinates
# A. Bus Stop data
bus_plot_data <- bus_long %>%
filter(is_weekend == 0) %>%
mutate(
X = st_coordinates(geometry)[,1],
Y = st_coordinates(geometry)[,2]
) %>%
st_drop_geometry()
# B. Crime data
crime_plot_data <- crime_sf %>%
mutate(
X = st_coordinates(geometry)[,1],
Y = st_coordinates(geometry)[,2]
) %>%
st_drop_geometry()
# 2. Create Maps
# Left: Ridership Density
p_ridership <- ggplot() +
geom_sf(data = philly_boundary, fill = "#f5f5f5", color = "grey80") +
stat_density_2d(
data = bus_plot_data,
aes(x = X, y = Y, fill = ..level.., weight = Ridership),
geom = "polygon",
alpha = 0.75
) +
scale_fill_distiller(palette = "Blues", direction = 1, guide = "none") +
labs(
title = "Weekday Ridership Hotspots",
subtitle = "High Transit Activity Zones"
) +
mapTheme
# Right：Crime Density
p_crime_map <- ggplot() +
geom_sf(data = philly_boundary, fill = "#f5f5f5", color = "grey80") +
stat_density_2d(
data = crime_plot_data,
aes(x = X, y = Y, fill = ..level..),
geom = "polygon",
alpha = 0.4,
bins = 30,
adjust = 0.5
) +
scale_fill_distiller(palette = "Reds", direction = 1, guide = "none") +
labs(
title = "Crime Hotspots",
subtitle = "High Incident Zones"
) +
mapTheme
# 3. Combine Side-by-Side
combined_map <- p_ridership + p_crime_map +
plot_annotation(
title = "Spatial Mismatch Analysis: Eyes on the Street vs. Targets?",
subtitle = "Left: Where people are (Ridership) | Right: Where crimes happen",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 12, color = "grey40", hjust = 0.5)
)
)
combined_map
Transit activity is anchored in Center City and radiates outward along the major arterials (Broad St.and Market St.). However, crime exhibits a polycentric distribution. Beyond the city center, we observe significant high-crime clusters in North and West Philadelphia (Kensington/Allegheny).
This spatial mismatch suggests that ridership volume alone cannot explain crime risk. While the downtown core attracts crime due to sheer foot traffic (‘Targets’), the peripheral hotspots are likely driven by other environmental factors—such as socioeconomic disadvantage or the presence of crime generators (e.g., alcohol outlets)—rather than transit volume alone.
bar_plot_data <- final_data %>%
st_drop_geometry() %>%
group_by(is_weekend_factor) %>%
summarise(
Avg_Ridership = mean(Ridership, na.rm = TRUE),
Avg_Crime_count = mean(Crime_Daily_Rate, na.rm = TRUE)
) %>%
# Long
pivot_longer(
cols = c(Avg_Ridership, Avg_Crime_count),
names_to = "Metric",
values_to = "Value"
) %>%
mutate(
Metric_Label = case_when(
Metric == "Avg_Ridership" ~ "Average Ridership",
Metric == "Avg_Crime_count" ~ "Average Crime Count"
)
)
ggplot(bar_plot_data, aes(x = is_weekend_factor, y = Value, fill = is_weekend_factor)) +
geom_col(width = 0.6, alpha = 0.9) +
facet_wrap(~Metric_Label, scales = "free_y") +
scale_fill_manual(values = c("Weekday" = "#3182bd", "Weekend" = "#de2d26")) +
labs(
title = "Volume vs. Risk: Weekday vs. Weekend",
subtitle = "Comparison of Average Ridership and Crime per Stop",
x = "",
y = "Average Count",
fill = "Time Period"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 10, color = "grey40", hjust = 0.5),
strip.text = element_text(size = 12, face = "bold"),
axis.text.x = element_text(size = 11, face = "bold"),
legend.position = "none" #
)
# Crime vs. Ridership (Interaction Plot)
ggplot(final_data, aes(x = Log_Ridership, y = Crime_Daily_Rate, color = is_weekend_factor)) +
geom_point(alpha = 0.1, size = 1) +
geom_smooth(method = "glm", method.args = list(family = "quasipoisson"), se = TRUE) +
scale_color_manual(values = c("Weekday" = "#3182bd", "Weekend" = "#de2d26")) +
labs(
title = "Does Ridership impact Crime differently on Weekends?",
subtitle = "Interaction Effect (Normalized by Number of Days)",
x = "Log(Daily Ridership)",
y = "Average Daily Crime Count (per 400m)",
color = "Time Period"
) +
plotTheme
Crime Follows the Crowd (Bar Chart):
The left bar charts show that crime volume drops significantly on weekends, mirroring the drop in ridership shown in the right chart. This confirms that crime is likely driven by opportunity—fewer people on the street simply means fewer targets and fewer conflicts.
The Risk “Conversion Rate” is Constant (Scatter Plot):
The interaction plot reveals a consistent positive correlation between ridership and crime counts across both weekdays and weekends. While the slopes exhibit only minor differences, the methodological value of this split is significant. By comparing the same bus stops during different time periods (Weekday vs. Weekend), we inherently control for static environmental factors, such as poverty rates, street lighting, and proximity to alcohol outlets, which remain constant regardless of the day.
Consequently, this acts as a quasi-experimental control: since the physical environment is fixed, the persistent upward trend suggests that ridership itself is a direct driver of crime risk. This supports the ‘Targets’ hypothesis over the ‘Eyes on the Street’ theory, indicating that higher passenger volumes attract opportunistic crime regardless of the temporal context.
library(patchwork)
library(scales)
plot_data <- final_data %>%
mutate(
log_crime = log(Crime_Daily_Rate + 0.01), # +0.01 to avoid log(0)
log_income = log(Avg_Income + 1)
)
# --- 1. POI & Infrastructure (Built Environment) ---
# A. Alcohol Outlets
p1 <- ggplot(plot_data, aes(x = Alcohol_Count, y = log_crime)) +
geom_point(alpha = 0.1, color = "#6A1B9A") +
geom_smooth(method = "loess", color = "red", se = FALSE) +
labs(title = "Log(Crime) vs. Alcohol Outlets",
subtitle = "Check for diminishing returns",
x = "Count of Alcohol Outlets", y = "Log(Crime Count)") +
plotTheme
# B. Street Lights
p2 <- ggplot(plot_data, aes(x = Light_Count, y = log_crime)) +
geom_point(alpha = 0.1, color = "#6A1B9A") +
geom_smooth(method = "loess", color = "red", se = FALSE) +
labs(title = "Log(Crime) vs. Street Lights",
subtitle = "Is the relationship linear?",
x = "Count of Street Lights", y = "") +
plotTheme
# C. Distance to Police
p3 <- ggplot(plot_data, aes(x = Dist_Police, y = log_crime)) +
geom_point(alpha = 0.1, color = "#6A1B9A") +
geom_smooth(method = "loess", color = "red", se = FALSE) +
labs(title = "Log(Crime) vs. Dist to Police",
subtitle = "Check for U-shape or threshold",
x = "Distance to Station (Miles)", y = "") +
plotTheme
# --- 2. Demographics (Census) ---
# D. Poverty Rate
p4 <- ggplot(plot_data, aes(x = Avg_Poverty, y = log_crime)) +
geom_point(alpha = 0.1, color = "#3182bd") +
geom_smooth(method = "loess", color = "orange", se = FALSE) +
scale_x_continuous(labels = scales::percent) +
labs(title = "Log(Crime) vs. Poverty Rate",
x = "Poverty Rate", y = "Log(Crime Count)") +
plotTheme
# E. Median Income
p5 <- ggplot(plot_data, aes(x = Avg_Income, y = log_crime)) +
geom_point(alpha = 0.1, color = "#3182bd") +
geom_smooth(method = "loess", color = "orange", se = FALSE) +
scale_x_continuous(labels = scales::dollar) +
labs(title = "Log(Crime) vs. Median Income",
subtitle = "Does wealth shield against crime?",
x = "Median Household Income", y = "") +
plotTheme
# F. Vacancy Rate
p6 <- ggplot(plot_data, aes(x = Avg_Vacancy, y = log_crime)) +
geom_point(alpha = 0.1, color = "#3182bd") +
geom_smooth(method = "loess", color = "orange", se = FALSE) +
scale_x_continuous(labels = scales::percent) +
labs(title = "Log(Crime) vs. Vacancy Rate",
x = "Housing Vacancy Rate", y = "") +
plotTheme
# --- Combine Plots ---
(p1 | p2 | p3) / (p4 | p5 | p6) +
plot_annotation(
title = "Exploratory Analysis: Variable Functional Forms",
subtitle = "Red Line = Loess Smoother (Non-linear trend)",
theme = theme(plot.title = element_text(size = 16, face = "bold"))
)
Before finalizing our model, we conducted an extensive exploratory analysis, testing relationships across multiple potential independent variables and their non-linear transformations. While we only present the variables with the most significant statistical relationships and policy implications here, this rigorous process of trial and error was essential to ensure the robustness of our final model. The charts above represent the distilled ‘risk signals’ identified from this comprehensive screening.
Alcohol Outlets: The plot shows that crime goes up quickly with the first few alcohol stores. However, as the number of stores gets very high, the line flattens out (Diminishing marginal returns). We will try using a log transformation for this feature. This tells the model that the first few stores have a bigger impact than the later ones.
Street Lights: The red trend line is not straight. It goes up as lights increase (usually in busy areas), but then it curves. In later this study, we used a polynomial (2nd degree) term. This allows the model to draw a curved line instead of a straight one to fit the data better.
Distance to Police: The line looks mostly straight and goes down slightly, so we kept this variable linear.
Poverty Rate: As poverty goes up, crime generally goes up. The line is mostly straight. We also kept it as a standard linear variable.
Median Income: The line looks a bit like a “U” shape. Crime is high in low-income areas, drops in middle-income areas, and rises slightly or stays flat in higher-income areas. We added a squared term (^2). This helps the model understand that crime can be high at both low and high income levels, but low in the middle.
Vacancy Rate: The line looks like an upside-down “U”. Crime goes up as vacancy increases, but when vacancy gets very high (empty neighborhoods), crime actually drops. We also added a squared term (^2) for this.
library(ggcorrplot)
numeric_vars <- final_data %>%
st_drop_geometry() %>%
dplyr::select(
Crime_Daily_Rate,
Ridership,
Alcohol_Count,
Light_Count,
Dist_Police,
Avg_Poverty,
Avg_Income,
Avg_Unemployment,
Avg_Vacancy
)
corr_matrix <- cor(numeric_vars, use = "complete.obs")
ggcorrplot(
corr_matrix,
method = "square",
type = "lower",
lab = TRUE,
lab_size = 3,
colors = c("#6D9EC1", "white", "#E46726"),
title = "Correlation Matrix of Features",
ggtheme = theme_minimal()
)
The correlation matrix reveals a critical issue of multicollinearity among the socioeconomic predictors. Specifically, Average Poverty Rate (Avg_Poverty) and Median Household Income (Avg_Income) exhibit a strong negative correlation of -0.81. Additionally, Unemployment Rate (Avg_Unemployment) shows a strong negative correlation with Income (-0.70) and a moderately high positive correlation with Poverty (0.66).
These high coefficients suggest that Poverty, Income, and Unemployment are capturing overlapping dimensions of neighborhood socioeconomic status. While this flags a risk of multicollinearity, correlation coefficients alone are insufficient for removing variables. Therefore, we will proceed to calculate the Variance Inflation Factor (VIF) in the next step to rigorously quantify the severity of multicollinearity and determine the appropriate strategy for feature selection or transformation.
ggplot(final_data, aes(x = is_weekend_factor, y = Crime_Daily_Rate, fill = is_weekend_factor)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white") +
coord_cartesian(ylim = c(0, quantile(final_data$Crime_Daily_Rate, 0.95))) +
scale_fill_manual(values = c("Weekday" = "#3182bd", "Weekend" = "#de2d26")) +
labs(
title = "Daily Crime Risk: Weekday vs. Weekend",
subtitle = "Comparison of Average Daily Crime Counts per Stop",
x = "Time Period",
y = "Average Daily Crime Count (per 400m)",
caption = "Note: Values represent daily averages to account for fewer weekend days per year."
) +
plotTheme
Stability of Crime Volume: The means (indicated by white diamonds) and medians (solid black lines) for both weekdays and weekends are nearly identical, appearing at approximately the same daily average level (~0.6 incidents). This pattern suggests that the aggregate demand for public safety resources at these bus stops does not fluctuate significantly throughout the week.
Justification for Pseudo-Panel Modeling： While this chart shows crime counts are constant, we know that generally ridership volume typically drops on weekends. If the outcome (crime) remains the same while the input (ridership) decreases, it implies that the risk intensity per rider is likely higher on weekends. Constructing a pseudo-panel allows us to mathematically capture this changing elasticity, rather than averaging it out.
model_1 <- glm.nb(Crime_Total_Count ~
Log_Ridership +
offset(log(Exposure_Days)),
data = final_data)
summary(model_1)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership + offset(log(Exposure_Days)),
data = final_data, init.theta = 1.962937689, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.328484 0.018925 -70.20 <2e-16 ***
Log_Ridership 0.243197 0.004944 49.19 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(1.9629) family taken to be 1)
Null deviance: 14335 on 11065 degrees of freedom
Residual deviance: 11852 on 11064 degrees of freedom
AIC: 92538
Number of Fisher Scoring iterations: 1
Theta: 1.9629
Std. Err.: 0.0274
2 x log-likelihood: -92531.6570 model_2 <- glm.nb(Crime_Total_Count ~
Log_Ridership * is_weekend_factor +
offset(log(Exposure_Days)),
data = final_data)
summary(model_2)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
offset(log(Exposure_Days)), data = final_data, init.theta = 1.968491375,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.434242 0.027417 -52.311 < 2e-16
Log_Ridership 0.263192 0.006798 38.716 < 2e-16
is_weekend_factorWeekend 0.186359 0.038085 4.893 9.92e-07
Log_Ridership:is_weekend_factorWeekend -0.034408 0.010040 -3.427 0.00061
(Intercept) ***
Log_Ridership ***
is_weekend_factorWeekend ***
Log_Ridership:is_weekend_factorWeekend ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(1.9685) family taken to be 1)
Null deviance: 14371 on 11065 degrees of freedom
Residual deviance: 11849 on 11062 degrees of freedom
AIC: 92509
Number of Fisher Scoring iterations: 1
Theta: 1.9685
Std. Err.: 0.0275
2 x log-likelihood: -92498.9340 model_3 <- glm.nb(Crime_Total_Count ~
Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) +
poly(Light_Count, 2) +
Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Income + I(Avg_Income^2) +
Avg_Unemployment +
Dist_Police +
offset(log(Exposure_Days)),
data = final_data)
summary(model_3)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) + poly(Light_Count, 2) + Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) + Avg_Income + I(Avg_Income^2) +
Avg_Unemployment + Dist_Police + offset(log(Exposure_Days)),
data = final_data, init.theta = 3.809131636, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -8.603e-01 8.710e-02 -9.877 < 2e-16
Log_Ridership 1.102e-01 5.485e-03 20.086 < 2e-16
is_weekend_factorWeekend 4.330e-02 2.965e-02 1.460 0.14416
log(Alcohol_Count + 1) 2.836e-01 7.490e-03 37.869 < 2e-16
poly(Light_Count, 2)1 2.074e+01 9.830e-01 21.094 < 2e-16
poly(Light_Count, 2)2 -2.401e+00 6.208e-01 -3.867 0.00011
Avg_Poverty 1.310e+00 1.157e-01 11.325 < 2e-16
Avg_Vacancy 5.836e-03 3.003e-03 1.943 0.05197
I(Avg_Vacancy^2) -1.034e-04 8.434e-05 -1.226 0.22030
Avg_Income -1.162e-05 1.435e-06 -8.097 5.64e-16
I(Avg_Income^2) 5.343e-11 7.797e-12 6.852 7.27e-12
Avg_Unemployment -1.235e-02 1.879e-03 -6.575 4.87e-11
Dist_Police -1.986e-01 1.149e-02 -17.278 < 2e-16
Log_Ridership:is_weekend_factorWeekend -2.319e-02 7.718e-03 -3.004 0.00266
(Intercept) ***
Log_Ridership ***
is_weekend_factorWeekend
log(Alcohol_Count + 1) ***
poly(Light_Count, 2)1 ***
poly(Light_Count, 2)2 ***
Avg_Poverty ***
Avg_Vacancy .
I(Avg_Vacancy^2)
Avg_Income ***
I(Avg_Income^2) ***
Avg_Unemployment ***
Dist_Police ***
Log_Ridership:is_weekend_factorWeekend **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.8091) family taken to be 1)
Null deviance: 25397 on 11065 degrees of freedom
Residual deviance: 11556 on 11052 degrees of freedom
AIC: 85737
Number of Fisher Scoring iterations: 1
Theta: 3.8091
Std. Err.: 0.0607
2 x log-likelihood: -85706.8380 model_4 <- glm.nb(Crime_Total_Count ~
Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) +
poly(Light_Count, 2) +
Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Income + I(Avg_Income^2) +
Avg_Unemployment +
Dist_Police +
factor(PSA_ID) +
offset(log(Exposure_Days)),
data = final_data,
control = glm.control(maxit = 100))
summary(model_4)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) + poly(Light_Count, 2) + Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) + Avg_Income + I(Avg_Income^2) +
Avg_Unemployment + Dist_Police + factor(PSA_ID) + offset(log(Exposure_Days)),
data = final_data, control = glm.control(maxit = 100), init.theta = 4.730969316,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.058e-01 1.071e-01 -3.789 0.000151
Log_Ridership 9.525e-02 5.201e-03 18.314 < 2e-16
is_weekend_factorWeekend 2.432e-05 2.758e-02 0.001 0.999296
log(Alcohol_Count + 1) 2.423e-01 8.067e-03 30.036 < 2e-16
poly(Light_Count, 2)1 2.468e+01 1.203e+00 20.510 < 2e-16
poly(Light_Count, 2)2 -5.178e+00 7.333e-01 -7.062 1.64e-12
Avg_Poverty 8.036e-01 1.364e-01 5.891 3.83e-09
Avg_Vacancy 9.710e-03 4.093e-03 2.372 0.017672
I(Avg_Vacancy^2) -1.480e-04 9.976e-05 -1.484 0.137899
Avg_Income -1.750e-05 1.684e-06 -10.396 < 2e-16
I(Avg_Income^2) 8.873e-11 8.819e-12 10.061 < 2e-16
Avg_Unemployment -1.928e-02 2.160e-03 -8.924 < 2e-16
Dist_Police -2.755e-01 1.534e-02 -17.954 < 2e-16
factor(PSA_ID)012 3.774e-02 7.062e-02 0.534 0.593118
factor(PSA_ID)021 5.770e-01 5.888e-02 9.800 < 2e-16
factor(PSA_ID)022 4.204e-01 6.138e-02 6.848 7.49e-12
factor(PSA_ID)023 2.167e-01 5.230e-02 4.144 3.42e-05
factor(PSA_ID)031 -2.512e-01 5.587e-02 -4.496 6.92e-06
factor(PSA_ID)032 1.861e-01 5.190e-02 3.586 0.000336
factor(PSA_ID)033 -1.804e-01 4.656e-02 -3.874 0.000107
factor(PSA_ID)051 -1.158e-01 6.476e-02 -1.788 0.073805
factor(PSA_ID)052 -3.001e-01 7.173e-02 -4.184 2.87e-05
factor(PSA_ID)053 -3.824e-01 8.832e-02 -4.330 1.49e-05
factor(PSA_ID)071 3.128e-01 6.760e-02 4.627 3.71e-06
factor(PSA_ID)072 -2.708e-01 6.293e-02 -4.303 1.69e-05
factor(PSA_ID)073 -3.697e-01 6.517e-02 -5.672 1.41e-08
factor(PSA_ID)081 4.798e-01 6.006e-02 7.989 1.36e-15
factor(PSA_ID)082 -1.388e-01 5.780e-02 -2.402 0.016302
factor(PSA_ID)083 -4.606e-03 6.210e-02 -0.074 0.940871
factor(PSA_ID)091 1.973e-01 5.179e-02 3.811 0.000139
factor(PSA_ID)092 1.249e-01 5.629e-02 2.219 0.026498
factor(PSA_ID)093 -2.243e-02 5.251e-02 -0.427 0.669180
factor(PSA_ID)094 4.995e-01 5.413e-02 9.229 < 2e-16
factor(PSA_ID)095 -1.978e-01 5.810e-02 -3.404 0.000663
factor(PSA_ID)121 4.666e-01 6.270e-02 7.443 9.85e-14
factor(PSA_ID)122 -1.427e-01 1.038e-01 -1.375 0.169163
factor(PSA_ID)123 8.914e-02 6.188e-02 1.441 0.149665
factor(PSA_ID)124 2.368e-01 7.720e-02 3.067 0.002159
factor(PSA_ID)141 -3.889e-02 5.571e-02 -0.698 0.485112
factor(PSA_ID)142 -1.625e-01 6.301e-02 -2.579 0.009916
factor(PSA_ID)143 1.960e-01 5.740e-02 3.415 0.000637
factor(PSA_ID)144 -1.241e-02 6.515e-02 -0.190 0.848975
factor(PSA_ID)151 1.524e-01 5.272e-02 2.891 0.003845
factor(PSA_ID)152 4.628e-01 5.253e-02 8.810 < 2e-16
factor(PSA_ID)153 5.232e-01 5.387e-02 9.713 < 2e-16
factor(PSA_ID)161 -3.663e-02 5.725e-02 -0.640 0.522209
factor(PSA_ID)162 -5.668e-02 5.983e-02 -0.947 0.343485
factor(PSA_ID)171 -1.610e-01 5.428e-02 -2.965 0.003022
factor(PSA_ID)172 -5.680e-01 6.898e-02 -8.234 < 2e-16
factor(PSA_ID)173 -3.975e-01 5.909e-02 -6.727 1.73e-11
factor(PSA_ID)181 -1.636e-01 6.676e-02 -2.451 0.014265
factor(PSA_ID)182 3.990e-02 5.880e-02 0.679 0.497423
factor(PSA_ID)183 -1.393e-01 5.350e-02 -2.603 0.009229
factor(PSA_ID)191 1.965e-01 5.347e-02 3.674 0.000238
factor(PSA_ID)192 1.749e-02 5.921e-02 0.295 0.767637
factor(PSA_ID)193 2.983e-01 5.720e-02 5.215 1.84e-07
factor(PSA_ID)221 2.063e-01 5.666e-02 3.642 0.000271
factor(PSA_ID)222 -1.035e-02 6.510e-02 -0.159 0.873735
factor(PSA_ID)223 2.647e-01 5.970e-02 4.434 9.23e-06
factor(PSA_ID)224 1.685e-01 6.918e-02 2.435 0.014876
factor(PSA_ID)241 1.121e-01 6.299e-02 1.780 0.075073
factor(PSA_ID)242 3.243e-01 6.764e-02 4.795 1.63e-06
factor(PSA_ID)243 8.122e-02 6.635e-02 1.224 0.220954
factor(PSA_ID)251 1.819e-01 6.727e-02 2.704 0.006846
factor(PSA_ID)252 -2.377e-01 6.457e-02 -3.682 0.000232
factor(PSA_ID)253 -3.656e-01 6.789e-02 -5.385 7.25e-08
factor(PSA_ID)254 -5.605e-02 6.495e-02 -0.863 0.388107
factor(PSA_ID)261 -1.929e-01 6.498e-02 -2.969 0.002990
factor(PSA_ID)262 1.561e-01 5.853e-02 2.667 0.007655
factor(PSA_ID)263 -9.630e-02 5.027e-02 -1.916 0.055399
factor(PSA_ID)351 -3.793e-02 5.096e-02 -0.744 0.456679
factor(PSA_ID)352 6.887e-02 5.383e-02 1.279 0.200771
factor(PSA_ID)353 -2.201e-02 5.796e-02 -0.380 0.704186
factor(PSA_ID)391 1.119e-02 5.277e-02 0.212 0.832001
factor(PSA_ID)392 -1.532e-01 5.705e-02 -2.686 0.007233
factor(PSA_ID)393 1.888e-01 6.928e-02 2.725 0.006431
Log_Ridership:is_weekend_factorWeekend -1.720e-02 7.148e-03 -2.406 0.016114
(Intercept) ***
Log_Ridership ***
is_weekend_factorWeekend
log(Alcohol_Count + 1) ***
poly(Light_Count, 2)1 ***
poly(Light_Count, 2)2 ***
Avg_Poverty ***
Avg_Vacancy *
I(Avg_Vacancy^2)
Avg_Income ***
I(Avg_Income^2) ***
Avg_Unemployment ***
Dist_Police ***
factor(PSA_ID)012
factor(PSA_ID)021 ***
factor(PSA_ID)022 ***
factor(PSA_ID)023 ***
factor(PSA_ID)031 ***
factor(PSA_ID)032 ***
factor(PSA_ID)033 ***
factor(PSA_ID)051 .
factor(PSA_ID)052 ***
factor(PSA_ID)053 ***
factor(PSA_ID)071 ***
factor(PSA_ID)072 ***
factor(PSA_ID)073 ***
factor(PSA_ID)081 ***
factor(PSA_ID)082 *
factor(PSA_ID)083
factor(PSA_ID)091 ***
factor(PSA_ID)092 *
factor(PSA_ID)093
factor(PSA_ID)094 ***
factor(PSA_ID)095 ***
factor(PSA_ID)121 ***
factor(PSA_ID)122
factor(PSA_ID)123
factor(PSA_ID)124 **
factor(PSA_ID)141
factor(PSA_ID)142 **
factor(PSA_ID)143 ***
factor(PSA_ID)144
factor(PSA_ID)151 **
factor(PSA_ID)152 ***
factor(PSA_ID)153 ***
factor(PSA_ID)161
factor(PSA_ID)162
factor(PSA_ID)171 **
factor(PSA_ID)172 ***
factor(PSA_ID)173 ***
factor(PSA_ID)181 *
factor(PSA_ID)182
factor(PSA_ID)183 **
factor(PSA_ID)191 ***
factor(PSA_ID)192
factor(PSA_ID)193 ***
factor(PSA_ID)221 ***
factor(PSA_ID)222
factor(PSA_ID)223 ***
factor(PSA_ID)224 *
factor(PSA_ID)241 .
factor(PSA_ID)242 ***
factor(PSA_ID)243
factor(PSA_ID)251 **
factor(PSA_ID)252 ***
factor(PSA_ID)253 ***
factor(PSA_ID)254
factor(PSA_ID)261 **
factor(PSA_ID)262 **
factor(PSA_ID)263 .
factor(PSA_ID)351
factor(PSA_ID)352
factor(PSA_ID)353
factor(PSA_ID)391
factor(PSA_ID)392 **
factor(PSA_ID)393 **
Log_Ridership:is_weekend_factorWeekend *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(4.731) family taken to be 1)
Null deviance: 30304 on 11065 degrees of freedom
Residual deviance: 11532 on 10989 degrees of freedom
AIC: 83893
Number of Fisher Scoring iterations: 1
Theta: 4.7310
Std. Err.: 0.0792
2 x log-likelihood: -83736.6850 model_5 <- glm.nb(Crime_Total_Count ~
Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) +
poly(Light_Count, 2) +
Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Income + I(Avg_Income^2) +
Avg_Unemployment +
Dist_Police +
factor(PSA_ID) +
log(Crime_Daily_Rate_lag + 0.001) +
offset(log(Exposure_Days)),
data = final_data,
control = glm.control(maxit = 100))
summary(model_5)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) + poly(Light_Count, 2) + Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) + Avg_Income + I(Avg_Income^2) +
Avg_Unemployment + Dist_Police + factor(PSA_ID) + log(Crime_Daily_Rate_lag +
0.001) + offset(log(Exposure_Days)), data = final_data, control = glm.control(maxit = 100),
init.theta = 13.48398402, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.557e-01 7.508e-02 3.405 0.000661
Log_Ridership 3.097e-02 3.559e-03 8.701 < 2e-16
is_weekend_factorWeekend 7.513e-03 1.983e-02 0.379 0.704750
log(Alcohol_Count + 1) 6.111e-02 5.773e-03 10.587 < 2e-16
poly(Light_Count, 2)1 3.775e+00 8.603e-01 4.388 1.14e-05
poly(Light_Count, 2)2 -8.197e-01 4.967e-01 -1.650 0.098911
Avg_Poverty 1.727e-01 9.536e-02 1.811 0.070160
Avg_Vacancy -8.960e-03 2.840e-03 -3.155 0.001604
I(Avg_Vacancy^2) 2.615e-04 6.921e-05 3.779 0.000158
Avg_Income -8.633e-06 1.166e-06 -7.405 1.31e-13
I(Avg_Income^2) 4.994e-11 6.027e-12 8.286 < 2e-16
Avg_Unemployment -1.039e-02 1.539e-03 -6.750 1.48e-11
Dist_Police -9.557e-02 1.108e-02 -8.629 < 2e-16
factor(PSA_ID)012 -1.441e-01 5.067e-02 -2.844 0.004459
factor(PSA_ID)021 1.421e-01 3.993e-02 3.558 0.000373
factor(PSA_ID)022 6.487e-02 4.244e-02 1.529 0.126353
factor(PSA_ID)023 5.092e-02 3.636e-02 1.400 0.161366
factor(PSA_ID)031 -2.420e-02 3.758e-02 -0.644 0.519537
factor(PSA_ID)032 3.348e-02 3.519e-02 0.951 0.341369
factor(PSA_ID)033 -7.734e-02 3.133e-02 -2.469 0.013559
factor(PSA_ID)051 5.406e-02 4.560e-02 1.186 0.235803
factor(PSA_ID)052 -1.807e-01 5.194e-02 -3.479 0.000504
factor(PSA_ID)053 -3.275e-01 6.822e-02 -4.801 1.58e-06
factor(PSA_ID)071 2.336e-01 4.844e-02 4.822 1.42e-06
factor(PSA_ID)072 -2.252e-01 4.649e-02 -4.845 1.27e-06
factor(PSA_ID)073 -1.690e-01 4.871e-02 -3.469 0.000521
factor(PSA_ID)081 1.977e-01 4.246e-02 4.655 3.23e-06
factor(PSA_ID)082 -2.421e-01 4.216e-02 -5.742 9.34e-09
factor(PSA_ID)083 -9.608e-02 4.610e-02 -2.084 0.037172
factor(PSA_ID)091 -1.083e-01 3.486e-02 -3.108 0.001886
factor(PSA_ID)092 8.268e-02 3.771e-02 2.193 0.028343
factor(PSA_ID)093 3.994e-02 3.498e-02 1.142 0.253427
factor(PSA_ID)094 2.249e-01 3.584e-02 6.276 3.47e-10
factor(PSA_ID)095 -1.073e-01 3.890e-02 -2.758 0.005810
factor(PSA_ID)121 4.142e-01 4.334e-02 9.557 < 2e-16
factor(PSA_ID)122 -1.703e-01 7.205e-02 -2.363 0.018121
factor(PSA_ID)123 -1.704e-01 4.189e-02 -4.067 4.75e-05
factor(PSA_ID)124 2.026e-01 5.139e-02 3.943 8.05e-05
factor(PSA_ID)141 -1.294e-01 3.887e-02 -3.330 0.000868
factor(PSA_ID)142 -1.265e-01 4.355e-02 -2.905 0.003670
factor(PSA_ID)143 -1.477e-03 4.006e-02 -0.037 0.970583
factor(PSA_ID)144 -8.455e-02 4.711e-02 -1.795 0.072659
factor(PSA_ID)151 1.626e-02 3.617e-02 0.450 0.652963
factor(PSA_ID)152 1.268e-01 3.589e-02 3.532 0.000412
factor(PSA_ID)153 2.305e-01 3.677e-02 6.269 3.62e-10
factor(PSA_ID)161 1.339e-01 3.881e-02 3.449 0.000562
factor(PSA_ID)162 1.178e-01 4.045e-02 2.913 0.003581
factor(PSA_ID)171 -6.720e-02 3.659e-02 -1.836 0.066302
factor(PSA_ID)172 -4.248e-01 4.869e-02 -8.724 < 2e-16
factor(PSA_ID)173 -2.044e-01 3.989e-02 -5.124 2.99e-07
factor(PSA_ID)181 -1.544e-01 4.483e-02 -3.445 0.000572
factor(PSA_ID)182 1.001e-01 3.898e-02 2.569 0.010197
factor(PSA_ID)183 -8.624e-02 3.617e-02 -2.384 0.017118
factor(PSA_ID)191 2.030e-01 3.696e-02 5.493 3.94e-08
factor(PSA_ID)192 6.114e-02 3.944e-02 1.550 0.121118
factor(PSA_ID)193 2.517e-01 3.986e-02 6.315 2.71e-10
factor(PSA_ID)221 6.874e-02 3.797e-02 1.811 0.070214
factor(PSA_ID)222 7.326e-03 4.428e-02 0.165 0.868600
factor(PSA_ID)223 -6.076e-02 3.978e-02 -1.528 0.126626
factor(PSA_ID)224 -1.448e-01 4.634e-02 -3.124 0.001781
factor(PSA_ID)241 -4.814e-02 4.310e-02 -1.117 0.264063
factor(PSA_ID)242 -7.662e-03 4.520e-02 -0.169 0.865409
factor(PSA_ID)243 1.327e-01 4.562e-02 2.908 0.003640
factor(PSA_ID)251 1.919e-01 4.508e-02 4.256 2.08e-05
factor(PSA_ID)252 -4.972e-02 4.469e-02 -1.113 0.265880
factor(PSA_ID)253 -2.549e-01 4.660e-02 -5.469 4.52e-08
factor(PSA_ID)254 4.112e-02 4.372e-02 0.941 0.346951
factor(PSA_ID)261 -7.614e-02 4.393e-02 -1.733 0.083059
factor(PSA_ID)262 5.008e-02 3.908e-02 1.282 0.199990
factor(PSA_ID)263 -6.978e-02 3.389e-02 -2.059 0.039467
factor(PSA_ID)351 -9.185e-02 3.527e-02 -2.604 0.009201
factor(PSA_ID)352 7.575e-02 3.668e-02 2.065 0.038890
factor(PSA_ID)353 -1.149e-01 3.954e-02 -2.905 0.003667
factor(PSA_ID)391 -8.094e-02 3.676e-02 -2.202 0.027682
factor(PSA_ID)392 -7.300e-02 3.878e-02 -1.883 0.059766
factor(PSA_ID)393 4.386e-02 4.647e-02 0.944 0.345290
log(Crime_Daily_Rate_lag + 0.001) 7.084e-01 6.903e-03 102.622 < 2e-16
Log_Ridership:is_weekend_factorWeekend -8.622e-03 5.060e-03 -1.704 0.088411
(Intercept) ***
Log_Ridership ***
is_weekend_factorWeekend
log(Alcohol_Count + 1) ***
poly(Light_Count, 2)1 ***
poly(Light_Count, 2)2 .
Avg_Poverty .
Avg_Vacancy **
I(Avg_Vacancy^2) ***
Avg_Income ***
I(Avg_Income^2) ***
Avg_Unemployment ***
Dist_Police ***
factor(PSA_ID)012 **
factor(PSA_ID)021 ***
factor(PSA_ID)022
factor(PSA_ID)023
factor(PSA_ID)031
factor(PSA_ID)032
factor(PSA_ID)033 *
factor(PSA_ID)051
factor(PSA_ID)052 ***
factor(PSA_ID)053 ***
factor(PSA_ID)071 ***
factor(PSA_ID)072 ***
factor(PSA_ID)073 ***
factor(PSA_ID)081 ***
factor(PSA_ID)082 ***
factor(PSA_ID)083 *
factor(PSA_ID)091 **
factor(PSA_ID)092 *
factor(PSA_ID)093
factor(PSA_ID)094 ***
factor(PSA_ID)095 **
factor(PSA_ID)121 ***
factor(PSA_ID)122 *
factor(PSA_ID)123 ***
factor(PSA_ID)124 ***
factor(PSA_ID)141 ***
factor(PSA_ID)142 **
factor(PSA_ID)143
factor(PSA_ID)144 .
factor(PSA_ID)151
factor(PSA_ID)152 ***
factor(PSA_ID)153 ***
factor(PSA_ID)161 ***
factor(PSA_ID)162 **
factor(PSA_ID)171 .
factor(PSA_ID)172 ***
factor(PSA_ID)173 ***
factor(PSA_ID)181 ***
factor(PSA_ID)182 *
factor(PSA_ID)183 *
factor(PSA_ID)191 ***
factor(PSA_ID)192
factor(PSA_ID)193 ***
factor(PSA_ID)221 .
factor(PSA_ID)222
factor(PSA_ID)223
factor(PSA_ID)224 **
factor(PSA_ID)241
factor(PSA_ID)242
factor(PSA_ID)243 **
factor(PSA_ID)251 ***
factor(PSA_ID)252
factor(PSA_ID)253 ***
factor(PSA_ID)254
factor(PSA_ID)261 .
factor(PSA_ID)262
factor(PSA_ID)263 *
factor(PSA_ID)351 **
factor(PSA_ID)352 *
factor(PSA_ID)353 **
factor(PSA_ID)391 *
factor(PSA_ID)392 .
factor(PSA_ID)393
log(Crime_Daily_Rate_lag + 0.001) ***
Log_Ridership:is_weekend_factorWeekend .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(13.484) family taken to be 1)
Null deviance: 64740 on 11065 degrees of freedom
Residual deviance: 11871 on 10988 degrees of freedom
AIC: 76243
Number of Fisher Scoring iterations: 1
Theta: 13.484
Std. Err.: 0.307
2 x log-likelihood: -76085.254 vif(model_5) GVIF Df GVIF^(1/(2*Df))
Log_Ridership 2.140750 1 1.463130
is_weekend_factor 8.319162 1 2.884296
log(Alcohol_Count + 1) 3.845278 1 1.960938
poly(Light_Count, 2) 11.981538 2 1.860493
Avg_Poverty 9.159654 1 3.026492
Avg_Vacancy 23.286389 1 4.825597
I(Avg_Vacancy^2) 14.466404 1 3.803473
Avg_Income 98.553277 1 9.927400
I(Avg_Income^2) 63.899807 1 7.993735
Avg_Unemployment 3.897398 1 1.974183
Dist_Police 2.828785 1 1.681899
factor(PSA_ID) 2368.120093 63 1.063606
log(Crime_Daily_Rate_lag + 0.001) 2.737960 1 1.654678
Log_Ridership:is_weekend_factor 8.153167 1 2.855375Avg_Income) exhibited an extremely high Generalized VIF score of 9.92, indicating strong collinearity with Poverty Rate (Avg_Poverty). Including both variables introduces noise and destabilizes coefficient estimates, as they essentially measure the same underlying concept of economic disadvantage.Avg_Income and its quadratic term, retaining Avg_Poverty as the primary proxy for socioeconomic status due to its stronger theoretical link to crime risk in urban literature. Additionally, we simplified the Street Lights variable from a polynomial to a linear term (Light_Count). Given that the model includes robust spatial fixed effects (PSA_ID), the complex non-linear variation is already largely absorbed by the spatial controls, making a linear specification for infrastructure more parsimonious and interpretable.model_6 <- glm.nb(Crime_Total_Count ~
Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) +
Light_Count +
Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Unemployment +
Dist_Police +
factor(PSA_ID) +
log(Crime_Daily_Rate_lag + 0.001) +
offset(log(Exposure_Days)),
data = final_data,
control = glm.control(maxit = 100))
summary(model_6)
Call:
glm.nb(formula = Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) + Light_Count + Avg_Poverty + Avg_Vacancy +
I(Avg_Vacancy^2) + Avg_Unemployment + Dist_Police + factor(PSA_ID) +
log(Crime_Daily_Rate_lag + 0.001) + offset(log(Exposure_Days)),
data = final_data, control = glm.control(maxit = 100), init.theta = 13.34075513,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.939e-01 4.323e-02 -4.485 7.30e-06
Log_Ridership 3.134e-02 3.564e-03 8.793 < 2e-16
is_weekend_factorWeekend 8.872e-03 1.987e-02 0.447 0.655203
log(Alcohol_Count + 1) 6.276e-02 5.758e-03 10.899 < 2e-16
Light_Count 1.154e-04 3.655e-05 3.158 0.001590
Avg_Poverty 3.399e-01 7.377e-02 4.607 4.08e-06
Avg_Vacancy -6.985e-03 2.840e-03 -2.460 0.013899
I(Avg_Vacancy^2) 2.459e-04 6.941e-05 3.543 0.000396
Avg_Unemployment -9.685e-03 1.477e-03 -6.555 5.55e-11
Dist_Police -9.825e-02 1.095e-02 -8.971 < 2e-16
factor(PSA_ID)012 -1.422e-01 5.058e-02 -2.811 0.004946
factor(PSA_ID)021 1.509e-01 3.998e-02 3.773 0.000161
factor(PSA_ID)022 1.014e-01 4.228e-02 2.399 0.016438
factor(PSA_ID)023 6.396e-02 3.636e-02 1.759 0.078605
factor(PSA_ID)031 2.073e-02 3.584e-02 0.578 0.563028
factor(PSA_ID)032 5.317e-02 3.469e-02 1.533 0.125309
factor(PSA_ID)033 -8.258e-02 3.112e-02 -2.654 0.007952
factor(PSA_ID)051 6.217e-02 4.519e-02 1.376 0.168953
factor(PSA_ID)052 -1.710e-01 5.178e-02 -3.302 0.000961
factor(PSA_ID)053 -3.309e-01 6.830e-02 -4.845 1.27e-06
factor(PSA_ID)071 2.547e-01 4.842e-02 5.259 1.45e-07
factor(PSA_ID)072 -2.214e-01 4.646e-02 -4.764 1.89e-06
factor(PSA_ID)073 -1.700e-01 4.882e-02 -3.482 0.000498
factor(PSA_ID)081 2.072e-01 4.248e-02 4.877 1.08e-06
factor(PSA_ID)082 -2.493e-01 4.219e-02 -5.909 3.45e-09
factor(PSA_ID)083 -1.196e-01 4.575e-02 -2.615 0.008924
factor(PSA_ID)091 -9.835e-02 3.360e-02 -2.927 0.003423
factor(PSA_ID)092 1.325e-01 3.576e-02 3.705 0.000212
factor(PSA_ID)093 3.605e-02 3.468e-02 1.040 0.298572
factor(PSA_ID)094 2.110e-01 3.357e-02 6.287 3.25e-10
factor(PSA_ID)095 -3.208e-02 3.647e-02 -0.880 0.379119
factor(PSA_ID)121 4.456e-01 4.336e-02 10.275 < 2e-16
factor(PSA_ID)122 -1.271e-01 7.201e-02 -1.765 0.077630
factor(PSA_ID)123 -1.146e-01 4.141e-02 -2.767 0.005659
factor(PSA_ID)124 2.308e-01 5.113e-02 4.515 6.34e-06
factor(PSA_ID)141 -1.059e-01 3.856e-02 -2.747 0.006022
factor(PSA_ID)142 -8.819e-02 4.343e-02 -2.030 0.042325
factor(PSA_ID)143 2.082e-02 4.008e-02 0.519 0.603501
factor(PSA_ID)144 -5.119e-02 4.695e-02 -1.090 0.275505
factor(PSA_ID)151 4.567e-02 3.592e-02 1.271 0.203590
factor(PSA_ID)152 1.576e-01 3.578e-02 4.404 1.06e-05
factor(PSA_ID)153 2.389e-01 3.681e-02 6.490 8.58e-11
factor(PSA_ID)161 1.985e-01 3.806e-02 5.215 1.84e-07
factor(PSA_ID)162 1.781e-01 3.975e-02 4.480 7.48e-06
factor(PSA_ID)171 -4.983e-02 3.509e-02 -1.420 0.155534
factor(PSA_ID)172 -4.390e-01 4.839e-02 -9.073 < 2e-16
factor(PSA_ID)173 -2.201e-01 3.978e-02 -5.533 3.15e-08
factor(PSA_ID)181 -9.756e-02 4.423e-02 -2.206 0.027403
factor(PSA_ID)182 1.344e-01 3.864e-02 3.479 0.000503
factor(PSA_ID)183 -4.389e-02 3.495e-02 -1.256 0.209178
factor(PSA_ID)191 1.931e-01 3.687e-02 5.237 1.63e-07
factor(PSA_ID)192 1.307e-01 3.848e-02 3.396 0.000684
factor(PSA_ID)193 2.633e-01 3.985e-02 6.608 3.88e-11
factor(PSA_ID)221 1.218e-01 3.758e-02 3.242 0.001189
factor(PSA_ID)222 7.118e-02 4.376e-02 1.627 0.103789
factor(PSA_ID)223 -1.989e-02 3.955e-02 -0.503 0.615088
factor(PSA_ID)224 -1.497e-01 4.627e-02 -3.235 0.001215
factor(PSA_ID)241 1.896e-03 4.282e-02 0.044 0.964684
factor(PSA_ID)242 4.314e-02 4.490e-02 0.961 0.336712
factor(PSA_ID)243 1.449e-01 4.519e-02 3.206 0.001345
factor(PSA_ID)251 2.562e-01 4.454e-02 5.752 8.84e-09
factor(PSA_ID)252 8.936e-03 4.432e-02 0.202 0.840212
factor(PSA_ID)253 -1.611e-01 4.549e-02 -3.540 0.000399
factor(PSA_ID)254 1.288e-01 4.258e-02 3.024 0.002497
factor(PSA_ID)261 -2.606e-02 4.353e-02 -0.599 0.549301
factor(PSA_ID)262 6.579e-02 3.819e-02 1.723 0.084969
factor(PSA_ID)263 -4.964e-02 3.330e-02 -1.491 0.135995
factor(PSA_ID)351 -7.299e-02 3.527e-02 -2.070 0.038487
factor(PSA_ID)352 1.256e-01 3.624e-02 3.466 0.000528
factor(PSA_ID)353 -5.349e-02 3.886e-02 -1.377 0.168638
factor(PSA_ID)391 -4.033e-02 3.657e-02 -1.103 0.270083
factor(PSA_ID)392 -1.015e-02 3.821e-02 -0.266 0.790478
factor(PSA_ID)393 1.247e-01 4.552e-02 2.739 0.006163
log(Crime_Daily_Rate_lag + 0.001) 7.118e-01 6.898e-03 103.194 < 2e-16
Log_Ridership:is_weekend_factorWeekend -8.741e-03 5.073e-03 -1.723 0.084891
(Intercept) ***
Log_Ridership ***
is_weekend_factorWeekend
log(Alcohol_Count + 1) ***
Light_Count **
Avg_Poverty ***
Avg_Vacancy *
I(Avg_Vacancy^2) ***
Avg_Unemployment ***
Dist_Police ***
factor(PSA_ID)012 **
factor(PSA_ID)021 ***
factor(PSA_ID)022 *
factor(PSA_ID)023 .
factor(PSA_ID)031
factor(PSA_ID)032
factor(PSA_ID)033 **
factor(PSA_ID)051
factor(PSA_ID)052 ***
factor(PSA_ID)053 ***
factor(PSA_ID)071 ***
factor(PSA_ID)072 ***
factor(PSA_ID)073 ***
factor(PSA_ID)081 ***
factor(PSA_ID)082 ***
factor(PSA_ID)083 **
factor(PSA_ID)091 **
factor(PSA_ID)092 ***
factor(PSA_ID)093
factor(PSA_ID)094 ***
factor(PSA_ID)095
factor(PSA_ID)121 ***
factor(PSA_ID)122 .
factor(PSA_ID)123 **
factor(PSA_ID)124 ***
factor(PSA_ID)141 **
factor(PSA_ID)142 *
factor(PSA_ID)143
factor(PSA_ID)144
factor(PSA_ID)151
factor(PSA_ID)152 ***
factor(PSA_ID)153 ***
factor(PSA_ID)161 ***
factor(PSA_ID)162 ***
factor(PSA_ID)171
factor(PSA_ID)172 ***
factor(PSA_ID)173 ***
factor(PSA_ID)181 *
factor(PSA_ID)182 ***
factor(PSA_ID)183
factor(PSA_ID)191 ***
factor(PSA_ID)192 ***
factor(PSA_ID)193 ***
factor(PSA_ID)221 **
factor(PSA_ID)222
factor(PSA_ID)223
factor(PSA_ID)224 **
factor(PSA_ID)241
factor(PSA_ID)242
factor(PSA_ID)243 **
factor(PSA_ID)251 ***
factor(PSA_ID)252
factor(PSA_ID)253 ***
factor(PSA_ID)254 **
factor(PSA_ID)261
factor(PSA_ID)262 .
factor(PSA_ID)263
factor(PSA_ID)351 *
factor(PSA_ID)352 ***
factor(PSA_ID)353
factor(PSA_ID)391
factor(PSA_ID)392
factor(PSA_ID)393 **
log(Crime_Daily_Rate_lag + 0.001) ***
Log_Ridership:is_weekend_factorWeekend .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(13.3408) family taken to be 1)
Null deviance: 64299 on 11065 degrees of freedom
Residual deviance: 11874 on 10991 degrees of freedom
AIC: 76310
Number of Fisher Scoring iterations: 1
Theta: 13.341
Std. Err.: 0.303
2 x log-likelihood: -76157.840 vif(model_6) GVIF Df GVIF^(1/(2*Df))
Log_Ridership 2.133246 1 1.460563
is_weekend_factor 8.301509 1 2.881234
log(Alcohol_Count + 1) 3.798142 1 1.948882
Light_Count 5.541211 1 2.353978
Avg_Poverty 5.428573 1 2.329930
Avg_Vacancy 23.126903 1 4.809044
I(Avg_Vacancy^2) 14.459269 1 3.802535
Avg_Unemployment 3.561355 1 1.887155
Dist_Police 2.751623 1 1.658802
factor(PSA_ID) 467.133380 63 1.049992
log(Crime_Daily_Rate_lag + 0.001) 2.717994 1 1.648634
Log_Ridership:is_weekend_factor 8.139565 1 2.852992Log_Ridership, has a very low GVIF score of 2.13. This is well below the concern threshold (usually 5 or 10), which proves that our main finding that ridership drives crime, is reliable and not distorted by other factors.library(modelsummary)
library(ggplot2)
library(dplyr)
# 1. Modify Coef Map
coef_map_refined <- c(
"Log_Ridership" = "Ridership (Log)",
"is_weekend_factorWeekend" = "Weekend Effect",
"Log_Ridership:is_weekend_factorWeekend" = "Interaction: Ridership × Weekend",
"is_weekend_factorWeekend:Log_Ridership" = "Interaction: Ridership × Weekend",
"log(Alcohol_Count + 1)" = "Alcohol Outlets",
"Light_Count" = "Street Lights (Linear)",
"Avg_Poverty" = "Poverty Rate",
"Avg_Vacancy" = "Vacancy Rate",
"I(Avg_Vacancy^2)" = "Vacancy Rate (Sq)",
"Avg_Unemployment" = "Unemployment Rate",
"Dist_Police" = "Dist to Police Station",
"log(Crime_Daily_Rate_lag + 0.001)" = "Temporal Lag (Prev. Q)"
)
# 2. Mapping
p_comparison <- modelplot(
list(
"M1: Ridership" = model_1,
"M2: +Interact" = model_2,
"M3: +Env/Demo" = model_3,
"M4: +PSA Fixed" = model_4,
"M5: +Lag" = model_5,
"M6: Final" = model_6
),
coef_map = coef_map_refined,
coef_omit = "Intercept",
conf_level = 0.95,
size = 0.8
) +
geom_vline(xintercept = 0, color = "red", linetype = "dashed") +
labs(
title = "Model Evolution: Stability of Key Drivers",
subtitle = "Comparing coefficients across model iterations (M1-M6)",
x = "Effect Size (Coefficient Estimate)",
y = "",
caption = "Note: To facilitate direct comparison, this plot displays only the variables selected for the Final Refined Model (M6)."
) +
theme_minimal() +
scale_color_brewer(palette = "Dark2") +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14),
axis.text.y = element_text(size = 10, face = "bold", color = "black"),
panel.grid.minor = element_blank()
)
p_comparison
This plot visualizes our journey from a naive baseline to a robust final model. When reading this chart:
Model 1(Ridership Only): Without controlling for any other factors, ridership appears to be a massive driver of crime. However, this estimate is likely inflated (biased) because high-ridership stops are often located in dense, low-income areas. The model is mistakenly attributing the neighborhood effect solely to ridership.
Model 2(+Interaction): We introduce the Weekend Effect and the interaction term. The main Ridership coefficient increases slightly. This confirms that the relationship isn’t identical across the week. By separating weekends, we begin to see that the “baseline risk” shifts, validating the need for a pseudo-panel approach, though the core finding (ridership = risk) remains consistent due to its positive coefficient.
Model 3(+Env & Demo): A dramatic shift occurs. The Poverty Rate coefficient appears and is extremely high (purple dot far to the right), while the Ridership coefficient drops significantly (moving left). This is the “Reality Check.” Once we account for poverty, we realize that socioeconomic status is the primary driver of crime, not just the bus stop itself. Ridership is still a significant risk factor (the coefficient is still positive), but its impact is much smaller than Model 1 suggested.
Model 4(+PSA Fixed Effects): The Poverty Rate coefficient shrinks (moves left, pink dot) compared to Model 3. By adding PSA (Police Service Area) fixed effects, we control for unobserved neighborhood characteristics (like local culture or police presence). The model no longer relies solely on “Poverty” to explain crime clusters, making the estimates for specific variables like Alcohol Outlets and Ridership more precise and trustworthy.
Model 5(+Lag): The Temporal Lag variable appears at the top with a very high positive coefficient. The strongest predictor of future crime is past crime. Adding this variable absorbs a huge amount of variance, creating a very strict test for our other variables.
Model 6(Final Refined): Ridership remains positive and significant (approx 0.1), proving that even after controlling for everything (poverty, history, location), more passengers still equal more targets. Alcohol shows a stable positive link. Dist to Police shows a negative coefficient (meaning crime is higher closer to stations), but this is probably because there are bias in original dataset where crimes nearby the police stations are more likely to be recorded.
f1 <- Crime_Total_Count ~ Log_Ridership + offset(log(Exposure_Days))
f2 <- Crime_Total_Count ~ Log_Ridership * is_weekend_factor + offset(log(Exposure_Days))
f3 <- Crime_Total_Count ~ Log_Ridership * is_weekend_factor + log(Alcohol_Count + 1) + poly(Light_Count, 2) + Avg_Poverty + Avg_Vacancy + I(Avg_Vacancy^2) + Avg_Income + I(Avg_Income^2) + Avg_Unemployment + Dist_Police + offset(log(Exposure_Days))
f4 <- Crime_Total_Count ~ Log_Ridership * is_weekend_factor + log(Alcohol_Count + 1) + poly(Light_Count, 2) + Avg_Poverty + Avg_Vacancy + I(Avg_Vacancy^2) + Avg_Income + I(Avg_Income^2) + Avg_Unemployment + Dist_Police + factor(PSA_ID) + offset(log(Exposure_Days))
f5 <- Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) + poly(Light_Count, 2) +
Avg_Poverty + Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Income + I(Avg_Income^2) + Avg_Unemployment +
Dist_Police + factor(PSA_ID) +
log(Crime_Daily_Rate_lag + 0.001) +
offset(log(Exposure_Days))
f6 <- Crime_Total_Count ~ Log_Ridership * is_weekend_factor +
log(Alcohol_Count + 1) +
Light_Count +
Avg_Poverty +
Avg_Vacancy + I(Avg_Vacancy^2) +
Avg_Unemployment +
Dist_Police +
factor(PSA_ID) +
log(Crime_Daily_Rate_lag + 0.001) +
offset(log(Exposure_Days))
# 2. Train 6 models
m1 <- glm.nb(f1, data = final_data)
m2 <- glm.nb(f2, data = final_data)
m3 <- glm.nb(f3, data = final_data)
m4 <- glm.nb(f4, data = final_data)
m5 <- glm.nb(f5, data = final_data)
m6 <- glm.nb(f6, data = final_data)
# 3. Gather Predictions
plot_data <- final_data %>%
dplyr::select(Crime_Total_Count) %>%
mutate(
Model_1_Pred = predict(m1, type = "response"),
Model_2_Pred = predict(m2, type = "response"),
Model_3_Pred = predict(m3, type = "response"),
Model_4_Pred = predict(m4, type = "response"),
Model_5_Pred = predict(m5, type = "response"),
Model_6_Pred = predict(m6, type = "response"),
) %>%
pivot_longer(
cols = starts_with("Model"),
names_to = "Model",
values_to = "Predicted_Count"
)
# 4. Map：Facet Wrap compare 6 models
ggplot(plot_data, aes(x = Predicted_Count, y = Crime_Total_Count)) +
geom_point(alpha = 0.2, color = "#3182bd", size = 0.8) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
facet_wrap(~Model, ncol = 3) +
labs(
title = "Predicted vs. Actual Crime Counts",
subtitle = "Comparing Model Fit: Final Model (M6) shows robust predictions",
x = "Predicted Crime Count",
y = "Actual Crime Count"
) +
theme_bw() +
coord_cartesian(xlim = c(0, 50), ylim = c(0, 50))
This plot also visualizes our journey from a naive baseline to a robust final model.
Overall, as more contextual and more spatial information are added, the predicted values align more closely with the 45° reference line, indicating improved model fit.
Models 1–2 rely mainly on ridership and weekend interaction, resulting in wide dispersion and systematic under-prediction at higher crime levels.
Models 3–4 incorporate environmental (alcohol outlets, lighting) and demographic features, producing a noticeably tighter prediction band and capturing more variation across stops.
Model 5 adds temporal information (lagged crime rate), which further stabilizes predictions and reduces bias.
Compared with Model 5, Model 6 trims several weak or redundant predictors while keeping the key effects. This parsimonious specification makes the model more stable and interpretable and helps reduce potential multicollinearity among highly correlated socioeconomic variables.
library(caret)
# 1. Set 5-Fold Cross Validation
set.seed(999)
folds <- createFolds(final_data$Crime_Total_Count, k = 5, list = TRUE)
run_cv <- function(formula, data, folds) {
mae_list <- c()
rmse_list <- c()
for (i in 1:length(folds)) {
# Split Data
test_idx <- folds[[i]]
train_set <- data[-test_idx, ]
test_set <- data[test_idx, ]
# Train Model
model <- tryCatch({
glm.nb(formula, data = train_set)
}, error = function(e) return(NULL))
if(!is.null(model)) {
preds <- predict(model, newdata = test_set, type = "response")
# Calculate Errors
actuals <- test_set$Crime_Total_Count
mae_list <- c(mae_list, mean(abs(actuals - preds)))
rmse_list <- c(rmse_list, sqrt(mean((actuals - preds)^2)))
}
}
return(c(mean(mae_list), mean(rmse_list)))
}
# 2. CV 6 models
results_m1 <- run_cv(f1, final_data, folds)
results_m2 <- run_cv(f2, final_data, folds)
results_m3 <- run_cv(f3, final_data, folds)
results_m4 <- run_cv(f4, final_data, folds)
results_m5 <- run_cv(f5, final_data, folds)
results_m6 <- run_cv(f6, final_data, folds)
# 3. Generate chart
validation_summary <- data.frame(
Model = c("1. Ridership Only", "2. +Interaction", "3. +Env & Demo", "4. +Spatial Fixed", "5. +Temporal lag", "6. Refined"),
MAE = c(results_m1[1], results_m2[1], results_m3[1], results_m4[1], results_m5[1], results_m6[1]),
RMSE = c(results_m1[2], results_m2[2], results_m3[2], results_m4[2], results_m5[2], results_m6[2])
) %>%
mutate(
Improvement_MAE = (MAE[1] - MAE) / MAE[1]
)
kable(validation_summary, digits = 3, caption = "5-Fold Cross-Validation Metrics") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE) %>%
column_spec(4, color = "green", bold = TRUE)| Model | MAE | RMSE | Improvement_MAE |
|---|---|---|---|
| 1. Ridership Only | 17.268 | 28.877 | 0.000 |
| 2. +Interaction | 17.209 | 28.868 | 0.003 |
| 3. +Env & Demo | 12.777 | 21.039 | 0.260 |
| 4. +Spatial Fixed | 11.409 | 18.343 | 0.339 |
| 5. +Temporal lag | 7.453 | 11.629 | 0.568 |
| 6. Refined | 7.475 | 11.662 | 0.567 |
These findings are validated through 5-fold cross-validation, which evaluates each model on unseen data rather than relying solely on in-sample fit. The consistent reduction in MAE and RMSE across folds confirms that the improvements are not artifacts of overfitting but reflect real gains in out-of-sample predictive performance.
Model 1(Ridership Only)
The model using only ridership intensity establishes a basic relationship between human activity levels and crime, but its predictive capacity is limited. Ridership alone captures general exposure but cannot explain spatial heterogeneity across stops.
Model 2(+ Interaction)
Adding the weekday–weekend interaction improves the model slightly.
Model 3(+ Environment&Demographic)
Alcohol outlet density, Lighting conditions poverty, income, vacancy, and unemployment: these variables substantially reduce prediction error, indicating that crime near transit stops is strongly shaped by both environmental risk factors and neighborhood disadvantage.
Model 4(+Spatial Fixed)
Introducing police station distance and especially PSA fixed effects further improves performance.
This suggests that:Crime patterns are partly structured by local policing jurisdictions and there are unobserved spatial factors (culture, enforcement style, land use patterns) that are constant within each PSA.PSA fixed effects absorb these stable spatial characteristics, making the model more robust.
Model 5(Temporal Lag)
Including lagged crime rate (previous quarter) yields the largest single improvement. Lagged crime captures temporal persistence—areas that experienced more crime in the past tend to remain active in the present.This is a very strong and stable predictor, dramatically improving accuracy.
Model 6(Refined)
Model 6 streamlines the specification by removing weaker or collinear variables.Despite using fewer predictors, its accuracy remains nearly identical to Model 5.
This indicates that a refined, parsimonious model with reduced multicollinearity can maintain strong predictive performance while improving interpretability and stability.
best_model <- model_6
final_data$resid_pearson <- residuals(best_model, type = "pearson")
final_data$resid_deviance <- residuals(best_model, type = "deviance")
# Spatial Weights Matrix
coords <- st_coordinates(final_data)
neighbor_nb <- knn2nb(knearneigh(coords, k = 8))
spatial_weights <- nb2listw(neighbor_nb, style = "W")
# Moran's I test
moran_result <- moran.test(final_data$resid_pearson, spatial_weights)
print(moran_result)
Moran I test under randomisation
data: final_data$resid_pearson
weights: spatial_weights
Moran I statistic standard deviate = 100.73, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
4.518144e-01 -9.037506e-05 2.012726e-05 Although we added spatial fixed effect(which is also a way to reduce spatial autocorrelation), the Moran’s I test on the model residuals shows a strong and highly significant level of spatial autocorrelation (Moran’s I = 0.452, p < 2.2 × 10⁻¹⁶).
This indicates that the residuals are not randomly distributed across space and that the model likely fails to capture important spatial dependence structures in the data.
In practical terms, nearby transit stops tend to have systematically similar over or under predictions, suggesting the presence of spatially clustered unobserved factors such as local policing practices, neighborhood context, or land-use patterns that are not fully absorbed even after including PSA fixed effects and temporal lag terms.
Given the strength of the residual spatial autocorrelation, the results imply that the model could benefit from incorporating explicit spatial components, such as:
Overall, the Moran’s I result shows that while the model fits well in non-spatial dimensions, unmodeled spatial processes remain significant and should be considered in future extensions.
# Spatial Distribution of Residuals
final_data$spatial_resid <- residuals(best_model, type = "deviance")
ggplot() +
geom_sf(data = philly_boundary, fill = "grey95", color = NA) +
geom_sf(data = final_data,
aes(color = spatial_resid),
size = 0.8, alpha = 0.7) +
scale_color_gradient2(
low = "blue", mid = "grey90", high = "red",
midpoint = 0,
name = "Deviance\nResidual",
limits = c(-3, 3),
oob = scales::squish
) +
labs(
title = "Map of Model Residuals",
subtitle = "Red = Unexpectedly High Crime (Under-predicted)\nBlue = Unexpectedly Low Crime (Over-predicted)"
) +
mapTheme
Red spots (positive deviance residuals):These locations experienced more crime than the model expected.
Blue spots (negative deviance residuals):These locations show less crime than predicted.
Although the residuals exhibit noticeable spatial clustering, which was confirmed by the significant Moran’s I test, the pattern does not form a clear, interpretable spatial structure.
The clusters of over and under prediction appear scattered across different parts of the city without aligning with obvious geographic boundaries, transit corridors, or neighborhood divisions.
This suggests that while some localized spatial dependence remains in the model, these patterns do not translate into a single coherent spatial process.
In other words, the remaining spatial signal is weakly structured and heterogeneous, making it difficult to summarize into a simple spatial rule or trend.
best_model <- model_6
model_data <- data.frame(
Fitted = fitted(best_model),
Residuals = residuals(best_model, type = "pearson")
)
p_resid_fitted <- ggplot(model_data, aes(x = log(Fitted), y = Residuals)) +
geom_point(alpha = 0.3, color = "#6A1B9A", size = 1.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
geom_smooth(method = "loess", color = "black", se = FALSE, linewidth = 0.8) +
labs(
title = "Residuals vs Fitted Values",
subtitle = "Checking for systematic bias in the Count Model",
x = "Log(Fitted Values) - Predicted Crime Count",
y = "Pearson Residuals"
) +
plotTheme
p_resid_fitted
This plot shows the residuals against the predicted crime counts from the final negative binomial model. Overall, the residuals are centered around zero without clear curvature or funnel shapes, suggesting that the model does not suffer from major misspecification or heteroskedasticity.
The “striped” patterns come from the integer nature of crime counts which is an effect commonly called the discreteness pattern in count models. This is not a model problem; it simply reflects that many stations share the same small crime counts , leading to stacked fitted values. The smooth LOESS line stays close to zero across the range, indicating no strong systematic bias.
In short, the residual plot shows that the negative binomial model is behaving as expected for count data, and there is no major evidence of nonlinearity or missing predictors.
p_qq <- ggplot(model_data, aes(sample = Residuals)) +
stat_qq(color = "#6A1B9A", size = 1.5, alpha = 0.5) +
stat_qq_line(color = "red", linetype = "dashed", linewidth = 1) +
labs(
title = "Q-Q Plot of Pearson Residuals",
subtitle = "Checking for extreme outliers in count data",
x = "Theoretical Quantiles",
y = "Sample Quantiles"
) +
plotTheme
p_qq
The Q–Q plot shows that the Pearson residuals follow the theoretical quantiles closely in the middle range, meaning the model captures the main distribution of crime counts well. The deviation at the upper tail indicates a few unusually high-crime stops that the model cannot fully explain—this is common in count data, especially when rare but extreme events occur.
Negative Binomial models are designed to handle overdispersion, and the relatively small tail deviation suggests that any remaining overdispersion is limited to a small number of extreme observations rather than a systematic model failure.
Overall, the Q–Q plot supports that the model fits the bulk of the data well, with only mild deviations at the extremes. Binomial assumption might still be struggling with extreme overdispersion.
library(dplyr)
library(ggplot2)
library(knitr)
library(kableExtra)
model_final <- model_6
coef_name_base <- "Log_Ridership"
coef_name_interact <- "Log_Ridership:is_weekend_factorWeekend"
beta_base <- coef(model_final)[coef_name_base]
beta_interact <- coef(model_final)[coef_name_interact]
vcov_matrix <- vcov(model_final)
# Calculate real slope in weekends and weekdays
# A. Weekday: Base
slope_weekday <- beta_base
se_weekday <- sqrt(vcov_matrix[coef_name_base, coef_name_base])
p_weekday <- 2 * (1 - pnorm(abs(slope_weekday / se_weekday)))
# B. Weekend: Base + Interaction
slope_weekend <- beta_base + beta_interact
# Calculate Var: Var(A+B) = Var(A) + Var(B) + 2*Cov(A,B)
var_weekend <- vcov_matrix[coef_name_base, coef_name_base] +
vcov_matrix[coef_name_interact, coef_name_interact] +
2 * vcov_matrix[coef_name_base, coef_name_interact]
se_weekend <- sqrt(var_weekend)
p_weekend <- 2 * (1 - pnorm(abs(slope_weekend / se_weekend)))
hypothesis_test <- data.frame(
Scenario = c("Weekday (Commuters)", "Weekend (Non-Routine)"),
Impact_Slope = c(slope_weekday, slope_weekend),
Std_Error = c(se_weekday, se_weekend),
P_Value = c(p_weekday, p_weekend)
) %>%
mutate(
Interpretation = case_when(
Impact_Slope > 0 & P_Value < 0.05 ~ "Positive & Significant (Target Hypothesis Supported)",
Impact_Slope < 0 & P_Value < 0.05 ~ "Negative & Significant (Eyes on Street Supported)",
TRUE ~ "Not Significant"
),
Impact_Slope = round(Impact_Slope, 3),
Std_Error = round(Std_Error, 3),
P_Value = ifelse(P_Value < 0.001, "< 0.001", round(P_Value, 3))
)
kbl(hypothesis_test, caption = "The Verdict: Does Ridership Drive Crime in Both Contexts?") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>%
row_spec(1:2, bold = TRUE, color = "black", background = "#e6f5ff") %>%
footnote(general = "Slopes > 0 indicate that higher ridership is associated with higher crime counts.")| Scenario | Impact_Slope | Std_Error | P_Value | Interpretation |
|---|---|---|---|---|
| Weekday (Commuters) | 0.031 | 0.004 | < 0.001 | Positive & Significant (Target Hypothesis Supported) |
| Weekend (Non-Routine) | 0.023 | 0.004 | < 0.001 | Positive & Significant (Target Hypothesis Supported) |
| Note: | ||||
| Slopes > 0 indicate that higher ridership is associated with higher crime counts. |
## 7.2.1 Top 50 Under-Policed stops
final_data <- final_data %>%
mutate(
Predicted_Crime = predict(best_model, type = "response"),
Resid_Raw = Crime_Total_Count - Predicted_Crime
)
top_50_anomalies <- final_data %>%
arrange(desc(Resid_Raw)) %>%
slice(1:50)
ggplot() +
geom_sf(data = philly_boundary, fill = "grey95", color = "white") +
geom_sf(data = final_data, color = "grey80", size = 0.5, alpha = 0.3) +
geom_sf(data = top_50_anomalies,
aes(size = Resid_Raw),
color = "steelblue",
alpha = 0.8) +
scale_size_continuous(name = "Excess Crimes\n(Actual - Predicted)") +
labs(
title = "Top 50 Under-policed Stops",
subtitle = "Locations where actual crime significantly exceeds model predictions.",
caption = "Blue dots represent stops performing worse than their environment, suggesting there are more police in need."
) +
mapTheme +
theme(legend.position = "right")
## 7.2.2 Top 50 Over-Policed stops
final_data <- final_data %>%
mutate(
Predicted_Crime = predict(best_model, type = "response"),
Resid_Raw = Predicted_Crime - Crime_Total_Count
)
top_50_anomalies <- final_data %>%
arrange(desc(Resid_Raw)) %>%
slice(1:50)
ggplot() +
geom_sf(data = philly_boundary, fill = "grey95", color = "white") +
geom_sf(data = final_data, color = "grey80", size = 0.5, alpha = 0.3) +
geom_sf(data = top_50_anomalies,
aes(size = Resid_Raw),
color = "#d73027",
alpha = 0.8) +
scale_size_continuous(name = "Excess Crimes\n(Predicted - Actual)") +
labs(
title = "Top 50 'Over-policed' Stops",
subtitle = "Locations where actual crime significantly under model predictions.",
caption = "Red dots represent stations which are statistically risky, but the actual crime counts are actually small. "
) +
mapTheme +
theme(legend.position = "right")
Instead of simply deploying police to the busiest stations (which is inefficient) or the poorest neighborhoods (which reinforces bias), our model offers a “Precision Policing” approach based on spatial anomalies. By calculating the residuals—the difference between actual and predicted crime—we spatially distinguished areas of resource saturation from areas of critical need.
As shown in our anomaly maps, the “Over-policed” stops (Red Dots) are heavily clustered in Center City, indicating that current security measures there are suppressing crime below predicted levels. In sharp contrast, we identified the “Top 50 Under-policed Stops” (Blue Dots). These blue anomalies represent locations where actual crime exceeds our model’s predictions by up to 90 incidents—a gap that local environmental factors like poverty or lighting cannot explain. We recommend SEPTA Transit Police utilize this “Hit List” to reallocate patrol units from the saturated downtown core to these specific outlying hotspots. Focusing resources on these specific anomalies yields the highest marginal return on public safety investment.
## 7.2.3 Top 10 High-Risk Anomaly Table
library(kableExtra)
top_10_table <- final_data %>%
mutate(
Predicted = predict(best_model, type = "response"),
Residual = Crime_Total_Count - Predicted
) %>%
arrange(desc(Residual)) %>%
slice(1:10) %>%
dplyr::select(
"Bus Stop Name" = Stop,
"Police District (PSA)" = PSA_ID,
"Actual Crime" = Crime_Total_Count,
"Model Predicted" = Predicted,
"Unexplained Excess" = Residual
)
kbl(top_10_table, digits = 1, caption = "The 'Hit List': Top 10 Stops with Highest Unexplained Crime Risk") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(5, bold = TRUE, color = "white", background = "#d73027") %>%
footnote(general = " 'Unexplained Excess' = Actual Crime minus Predicted Crime based on local environment. Positive values indicate specific local security failures.")| Bus Stop Name | Police District (PSA) | Actual Crime | Model Predicted | Unexplained Excess | geometry |
|---|---|---|---|---|---|
| Locust St & 17th St | 093 | 295 | 198.2 | 96.8 | POINT (2691904 234766) |
| Whitby Av & 53rd St | 124 | 138 | 52.3 | 85.7 | POINT (2675321 233206.7) |
| 17th St & Locust St | 093 | 276 | 190.9 | 85.1 | POINT (2691903 234716.1) |
| 54th St & Willows Av | 124 | 142 | 58.2 | 83.8 | POINT (2675499 232607.7) |
| 52nd St & Jefferson St | 193 | 176 | 94.5 | 81.5 | POINT (2675925 245535.5) |
| 54th St & Whitby Av - FS | 124 | 136 | 55.4 | 80.6 | POINT (2675117 232954.9) |
| 54th St & Whitby Av | 124 | 137 | 57.2 | 79.8 | POINT (2675141 233000) |
| Whitby Av & 53rd St - FS | 124 | 133 | 53.5 | 79.5 | POINT (2675326 233279) |
| Jefferson St & 52nd St - FS | 193 | 176 | 98.9 | 77.1 | POINT (2675866 245528.4) |
| 52nd St & Heston St | 193 | 176 | 98.9 | 77.1 | POINT (2675998 245602.9) |
| Note: | |||||
| 'Unexplained Excess' = Actual Crime minus Predicted Crime based on local environment. Positive values indicate specific local security failures. |
A major concern in crime modeling is the potential to stigmatize disadvantaged communities. We addressed this by incorporating Census Tract Fixed Effects (Phase 4) and controlling for poverty rates. This ensures our model compares bus stops within the same neighborhood context, rather than unfairly comparing a stop in North Philadelphia to one in Center City. However, we acknowledge that our final model still exhibits some spatial autocorrelation (Moran’s I \(\approx\) 0.29), indicating that unobserved spatial clusters remain.
While the Negative Binomial model effectively handles the over-dispersed count data, future iterations should explore Bayesian Spatial Models to fully resolve the complex spatial dependencies inherent in urban crime data.